This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the Lindgren, Rue, and Lindström (2011) stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation. Finally, Log-Gaussian Cox process models are fit using the pseudodata approach of Simpson et al. (2016).
Blangiardo and Cameletti (2015) section 6.1.
Toy dataset from Blangiardo and Cameletti (2015). A spatial process is observed via a random sampling plan that concentrates observations in the lower left of the unit square.
# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')
# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))
# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)
# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')
# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))
# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))
# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)
# Fit the model with INLA.
toy_fit <- inla(
y ~ -1 + intercept + f(spatial_field, model = spde),
data = inla.stack.data(stacks),
control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)
# Output posterior summaries.
toy_fit$summary.fixed
toy_fit$summary.hyperpar
# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']
# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005). The illustrations below fit a stationary LGCP model with no covariates, a Matèrn covariance function with \(\nu = 1\), and INLA’s default (flat?) priors for all parameters.
# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')
# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
N_QUADS <- 10
QUAD_SIZE <- 50
w_edge <- Frame(bei)$xrange[1]
e_edge <- Frame(bei)$xrange[2]
s_edge <- Frame(bei)$yrange[1]
n_edge <- Frame(bei)$yrange[2]
botleft <- cbind(
runif(N_QUADS, w_edge, e_edge - QUAD_SIZE),
runif(N_QUADS, s_edge, n_edge - QUAD_SIZE)
)
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
cbind(
botleft[r, 1] + c(0, 0, QUAD_SIZE, QUAD_SIZE),
botleft[r, 2] + c(0, QUAD_SIZE, QUAD_SIZE, 0)
)
)})
bei_win <- do.call(
union.owin,
apply(botleft, 1, function(x){return(
owin(x[1] + c(0, QUAD_SIZE), x[2] + c(0, QUAD_SIZE))
)})
)
bei_hole <- bei[complement.owin(bei_win, frame = Frame(bei))]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)
plot(bei_hole, main = 'Observed Region with Holes', pch = '.', cols = 'black')
plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)
The mesh should be fine in observed regions to accurately represent complex local structure but can be coarsened in unobserved regions where the model cannot infer as much detail. I am including margins to explore model behavior away from the observed regions.
# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25
MAX_EDGE_EXT <- 50
MARGIN <- 100
# Mesh covering the site.
bei_boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
bei_full_mesh <- inla.mesh.create(
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_full_spde <- inla.spde2.matern(mesh = bei_full_mesh)
plot(bei_full_mesh, asp = 1, main = 'Fine Mesh')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh')
# Mesh including a margin outside the site.
margin_mesh <- inla.mesh.2d(
loc = bei_full_mesh$loc[,1:2], # Include nodes from site.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_spde <- inla.spde2.matern(mesh = margin_mesh)
plot(margin_mesh, asp = 1, main = 'Fine Mesh with Coarse Margin')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Margin')
# Meshs with coarser resolution in quadrats.
quad_hole <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = bei_interior[[i]], grp = i - 1))
})
)
bei_hole_mesh0 <- inla.mesh.create(
boundary = list(bei_boundary, quad_hole),
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_hole0_spde <- inla.spde2.matern(mesh = bei_hole_mesh0)
plot(bei_hole_mesh0, asp = 1, main = 'Fine Mesh with Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Holes')
bei_hole_mesh <- inla.mesh.create(
loc = bei_hole_mesh0$loc[,1:2], # Include nodes from mesh with holes.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_hole_spde <- inla.spde2.matern(mesh = bei_hole_mesh)
plot(bei_hole_mesh, asp = 1, main = 'Fine Mesh with Coarse Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Holes')
# Meshs with finer resolution in quadrats.
quad_bnd <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = apply(bei_interior[[i]], 2, rev), grp = i - 1))
})
)
bei_samp_mesh0 <- inla.mesh.create(
boundary = quad_bnd,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_samp0_spde <- inla.spde2.matern(mesh = bei_samp_mesh0)
plot(bei_samp_mesh0, asp = 1, main = 'Fine Mesh in Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh in Quadrats')
bei_samp_mesh <- inla.mesh.create(
loc = bei_samp_mesh0$loc[,1:2], # Include nodes from mesh in quads.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_samp_spde <- inla.spde2.matern(mesh = bei_samp_mesh)
plot(bei_samp_mesh, asp = 1, main = 'Coarse Mesh with Fine Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Coarse Mesh with Fine Quadrats')
# Meshes with varying resolution in quadrats and a margin.
margin_hole <- inla.mesh.2d(
loc = bei_hole_mesh$loc[,1:2], # Include nodes from mesh with holes.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_hole_spde <- inla.spde2.matern(mesh = margin_hole)
plot(margin_hole, asp = 1, main = 'Coarse Holes and Margin')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Coarse Holes and Margin')
margin_samp <- inla.mesh.2d(
loc = bei_samp_mesh$loc[,1:2], # Include nodes from quads.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_samp_spde <- inla.spde2.matern(mesh = margin_samp)
plot(margin_samp, asp = 1, main = 'Fine Quadrats and Coarse Margin')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Quadrats and Coarse Margin')
# Identify which mesh nodes are in the oberserved region
# and create SPDE projectors for spatial mapping.
NPIX_X <- 400
NPIX_Y <- 200
obs_full <- rep(0, margin_mesh$n)
obs_full[inla.over_sp_mesh(as(Window(bei), 'SpatialPolygons'), margin_mesh, 'vertex')] <- 1
proj_margin_mesh <- inla.mesh.projector(margin_mesh, dims = c(NPIX_X, NPIX_Y))
obs_hole <- rep(0, margin_hole$n)
obs_hole[inla.over_sp_mesh(as(Window(bei_hole), 'SpatialPolygons'), margin_hole, 'vertex')] <- 1
proj_margin_hole <- inla.mesh.projector(margin_hole, dims = c(NPIX_X, NPIX_Y))
obs_samp <- rep(0, margin_samp$n)
obs_samp[inla.over_sp_mesh(as(Window(bei_samp), 'SpatialPolygons'), margin_samp, 'vertex')] <- 1
proj_margin_samp <- inla.mesh.projector(margin_samp, dims = c(NPIX_X, NPIX_Y))
Count the points in grid cells and fit a Poisson GLMM with the observed area in the cell as an exposure variable. The plots of the cell coutns are blank (white) where cells have no observed area.
NGRID_X <- 40
NGRID_Y <- 20
centers <- gridcenters(
dilation(bei_window_full, max(NGRID_X, NGRID_Y)),
NGRID_X, NGRID_Y)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_df))){
bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
bei$x < bei_df$x[r] + dx &
bei$y >= bei_df$y[r] - dy &
bei$y < bei_df$y[r] + dy)
bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}
)
user system elapsed
0.448 0.012 0.463
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_df$area > 0, bei_df$count, NA), nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')
# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(
margin_mesh,
as.matrix(bei_df[bei_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
full_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_spde$n.spde)
full_stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = margin_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))
# Set up stacks for prediction.
full_stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(full_stack_index), tag = 'pred_latent')
full_stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
full_stacks <- inla.stack(full_stack_est, full_stack_latent, full_stack_response)
# Fit the model with INLA.
system.time(
bei_full_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(full_stacks),
control.predictor = list(A = inla.stack.A(full_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
900.984 0.464 121.568
# Output posterior summaries.
bei_full_fit$summary.fixed
bei_full_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, bei_full_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, bei_full_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
bei_hole_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_hole_df))){
bei_hole_df$count[r] <- sum(bei_hole$x >= bei_hole_df$x[r] - dx &
bei_hole$x < bei_hole_df$x[r] + dx &
bei_hole$y >= bei_hole_df$y[r] - dy &
bei_hole$y < bei_hole_df$y[r] + dy)
bei_hole_df$area[r] <- area(Window(bei_hole)[owin(c(bei_hole_df$x[r] - dx, bei_hole_df$x[r] + dx), c(bei_hole_df$y[r] - dy, bei_hole_df$y[r] + dy))])
}
)
user system elapsed
1.208 0.004 1.210
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_hole_df$area > 0, bei_hole_df$count, NA), nrow = length(unique(bei_hole_df$x)))), unique(bei_hole_df$x), unique(bei_hole_df$y), unitname = 'meters'), ncolcours = range(bei_hole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
hole_A_est <- inla.spde.make.A(
margin_hole,
as.matrix(bei_hole_df[bei_hole_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
hole_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_hole_spde$n.spde)
hole_stack_est <- inla.stack(data = list(count = bei_hole_df$count[bei_hole_df$area > 0], larea = log(bei_hole_df$area[bei_hole_df$area > 0])), A = list(hole_A_est), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
hole_A_pred <- inla.spde.make.A(mesh = margin_hole, loc = as.matrix(bei_hole_df[,c('x', 'y')]))
# Set up stacks for prediction.
hole_stack_latent <- inla.stack(data = list(xi = NA), A = list(hole_A_pred), effects = list(hole_stack_index), tag = 'pred_latent')
hole_stack_response <- inla.stack(data = list(count = NA), A = list(hole_A_pred), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
hole_stacks <- inla.stack(hole_stack_est, hole_stack_latent, hole_stack_response)
# Fit the model with INLA.
system.time(
bei_hole_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_hole_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(hole_stacks),
control.predictor = list(A = inla.stack.A(hole_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
827.364 0.460 112.739
# Output posterior summaries.
bei_hole_fit$summary.fixed
bei_hole_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, bei_hole_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, bei_hole_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
bei_samp_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_samp_df))){
bei_samp_df$count[r] <- sum(bei_samp$x >= bei_samp_df$x[r] - dx &
bei_samp$x < bei_samp_df$x[r] + dx &
bei_samp$y >= bei_samp_df$y[r] - dy &
bei_samp$y < bei_samp_df$y[r] + dy)
bei_samp_df$area[r] <- area(Window(bei_samp)[owin(c(bei_samp_df$x[r] - dx, bei_samp_df$x[r] + dx), c(bei_samp_df$y[r] - dy, bei_samp_df$y[r] + dy))])
}
)
user system elapsed
0.972 0.000 0.974
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_samp_df$area > 0, bei_samp_df$count, NA), nrow = length(unique(bei_samp_df$x)))), unique(bei_samp_df$x), unique(bei_samp_df$y), unitname = 'meters'), ncolcours = range(bei_samp_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
samp_A_est <- inla.spde.make.A(
margin_samp,
as.matrix(bei_samp_df[bei_samp_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
samp_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_samp_spde$n.spde)
samp_stack_est <- inla.stack(data = list(count = bei_samp_df$count[bei_samp_df$area > 0], larea = log(bei_samp_df$area[bei_samp_df$area > 0])), A = list(samp_A_est), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
samp_A_pred <- inla.spde.make.A(mesh = margin_samp, loc = as.matrix(bei_samp_df[,c('x', 'y')]))
# Set up stacks for prediction.
samp_stack_latent <- inla.stack(data = list(xi = NA), A = list(samp_A_pred), effects = list(samp_stack_index), tag = 'pred_latent')
samp_stack_response <- inla.stack(data = list(count = NA), A = list(samp_A_pred), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
samp_stacks <- inla.stack(samp_stack_est, samp_stack_latent, samp_stack_response)
# Fit the model with INLA.
system.time(
bei_samp_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_samp_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(samp_stacks),
control.predictor = list(A = inla.stack.A(samp_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
267.940 0.292 37.552
# Output posterior summaries.
bei_samp_fit$summary.fixed
bei_samp_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, bei_samp_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, bei_samp_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.
The pseudodata are constructed as follows.
Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to their respective integration weights times the intensity at that their locations (so the integration weight is an exposure variable), observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to
\[\prod_{i=1}^{p+n} \eta_{i}^{y_{i}} \exp(-\alpha_{i} \eta_{i}).\]
(Is there a missing \(\alpha_{i}\)?)
Known “sampling effort” is accomodated by scaling the integration weights by the probability that a point would have been observed at that node (given the sampling plan), i.e. nodes in observed regions have the \(\alpha_{i}\) defined above and nodes in unobserved regions have \(\alpha_{i} = 0\). Weights can be scaled by values other than 0 or 1 to account for e.g. distance sampling with a (known) detection function or a (known) false negative rate for detection equipment. This adjustment assumes the point process of interest is observed through a known thinning process and then allows inference back to the intensity function of the unthinned process. More complicated detection processes are possible, e.g. Yuan et al. (2017) fit an LGCP to data obtianed through distance sampling while using splines to model the detection function.
full_pts <- cbind(bei$x, bei$y)
# Contruct the SPDE A matrix for nodes and points.
full_nV <- margin_mesh$n
full_nData <- dim(full_pts)[1]
full_LocationMatrix <- inla.mesh.project(margin_mesh, full_pts)$A
full_IntegrationMatrix <- sparseMatrix(i = 1:full_nV, j = 1:full_nV, x = rep(1, full_nV))
full_ObservationMatrix <- rbind(full_IntegrationMatrix, full_LocationMatrix)
# Get the integration weights.
full_IntegrationWeights <- diag(inla.mesh.fem(margin_mesh)$c0)
full_E_point_process <- c(obs_full * full_IntegrationWeights, rep(0, full_nData))
# Create the psuedodata.
full_fake_data <- c(rep(0, full_nV), rep(1, full_nData))
# Fit model to full site.
full_formula <- y ~ -1 + intercept + f(idx, model = margin_spde) # No covariates.
full_data <- list(y = full_fake_data, idx = 1:full_nV, intercept = rep(1, full_nV))
system.time(
result_full <- inla(
formula = full_formula,
data = full_data,
family = 'poisson',
control.predictor = list(A = full_ObservationMatrix),
E = full_E_point_process,
verbose = TRUE
)
)
user system elapsed
422.472 0.512 61.984
result_full$summary.fixed
result_full$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
hole_pts <- cbind(bei_hole$x, bei_hole$y)
# Contruct the SPDE A matrix for nodes and points.
hole_nV <- margin_hole$n
hole_nData <- dim(hole_pts)[1]
hole_LocationMatrix <- inla.mesh.project(margin_hole, hole_pts)$A
hole_IntegrationMatrix <- sparseMatrix(i = 1:hole_nV, j = 1:hole_nV, x = rep(1, hole_nV))
hole_ObservationMatrix <- rbind(hole_IntegrationMatrix, hole_LocationMatrix)
# Get the integration weights.
hole_IntegrationWeights <- diag(inla.mesh.fem(margin_hole)$c0)
hole_E_point_process <- c(obs_hole * hole_IntegrationWeights, rep(0, hole_nData))
# Create the psuedodata.
hole_fake_data <- c(rep(0, hole_nV), rep(1, hole_nData))
# Fit model to site with holes.
hole_formula <- y ~ -1 + intercept + f(idx, model = margin_hole_spde) # No covariates.
hole_data <- list(y = hole_fake_data, idx = 1:hole_nV, intercept = rep(1, hole_nV))
system.time(
result_hole <- inla(
formula = hole_formula,
data = hole_data,
family = 'poisson',
control.predictor = list(A = hole_ObservationMatrix),
E = hole_E_point_process,
verbose = TRUE
)
)
user system elapsed
359.856 0.440 52.879
result_hole$summary.fixed
result_hole$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
samp_pts <- cbind(bei_samp$x, bei_samp$y)
# Contruct the SPDE A matrix for nodes and points.
samp_nV <- margin_samp$n
samp_nData <- dim(samp_pts)[1]
samp_LocationMatrix <- inla.mesh.project(margin_samp, samp_pts)$A
samp_IntegrationMatrix <- sparseMatrix(i = 1:samp_nV, j = 1:samp_nV, x = rep(1, samp_nV))
samp_ObservationMatrix <- rbind(samp_IntegrationMatrix, samp_LocationMatrix)
# Get the integration weights.
samp_IntegrationWeights <- diag(inla.mesh.fem(margin_samp)$c0)
samp_E_point_process <- c(obs_samp * samp_IntegrationWeights, rep(0, samp_nData))
# Create the psuedodata.
samp_fake_data <- c(rep(0, samp_nV), rep(1, samp_nData))
# Fit model to quadrat-sampled site.
samp_formula <- y ~ -1 + intercept + f(idx, model = margin_samp_spde) # No covariates.
samp_data <- list(y = samp_fake_data, idx = 1:samp_nV, intercept = rep(1, samp_nV))
system.time(
result_samp <- inla(
formula = samp_formula,
data = samp_data,
family = 'poisson',
control.predictor = list(A = samp_ObservationMatrix),
E = samp_E_point_process,
verbose = TRUE
)
)
user system elapsed
296.832 0.648 45.128
result_samp$summary.fixed
result_samp$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
inlabruinlabru has a wrapper function to fit an LGCP with INLA, with a relatively easy-to-use interface for defining models and predicting arbitrary functions of latent variables. However, it is poorly documented, slow, and the documentation seems to imply that it does not support variable sampling effort (even though this appears to work).
bei_full_spdf <- as.SpatialPoints.ppp(bei)
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = margin_spde) + Intercept
system.time(
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf, E = obs_full, options = list(verbose = TRUE))
)
user system elapsed
3413.376 1.468 455.392
bei_full_lgcp$summary.fixed
bei_full_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_full <- predict(bei_full_lgcp, pixels(margin_mesh), ~ mySmooth)
plot(lambda_full, main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
plot(lambda_full['sd'] ,main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = margin_hole_spde) + Intercept
system.time(
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf, E = obs_hole, options = list(verbose = TRUE))
)
user system elapsed
3021.376 1.240 391.739
bei_hole_lgcp$summary.fixed
bei_hole_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_hole <- predict(bei_hole_lgcp, pixels(margin_hole), ~ mySmooth)
plot(lambda_hole, main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
plot(lambda_hole['sd'], main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = margin_samp_spde) + Intercept
system.time(
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf, E = obs_samp, options = list(verbose = TRUE))
)
user system elapsed
264.924 0.740 45.721
bei_samp_lgcp$summary.fixed
bei_samp_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_samp <- predict(bei_samp_lgcp, pixels(margin_samp), ~ mySmooth)
plot(lambda_samp, main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
plot(lambda_samp['sd'], main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
All three methods give similar results for the full dataset and the dataset with holes, even when the gridding method uses a coarse grid. The intercept and random effects are shifted from method to method, but these are not separately indentifiable and the shifts cancel each other out. The methods each have different artifacts and edge effects apparent in the results from the sampled dataset. The pseudodata approach is the fastest except when a very coarse grid is used and a small region is observed.
Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.
Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.
Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.
Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.
Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.
Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.
Yuan, Yuan, Fabian E Bachl, Finn Lindgren, David L Borchers, Janine B Illian, Stephen T Buckland, Haavard Rue, Tim Gerrodette, and others. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” The Annals of Applied Statistics 11 (4): 2270–97.